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Abstract 

We analyze the autocorrelations for the Local Hybrid Monte Carlo algorithm [1] in the 
context of free field theory. In this case this is just Adler's overrelaxation algorithm 
[2]. We consider the algorithm with even/odd, lexicographic, and random updates, and 
show that its efficiency depends crucially on this ordering of sites when optimized for a 
given class of operators. In particular, we show that, contrary to previous expectations, 
it is possible to eliminate critical slowing down (zim = 0) for a class of interesting 
observables, including the magnetic susceptibility: this can be done with lexicographic 
updates but is not possible with even/odd (zi^t = 1) or random (zi^t = 2) updates. 
We are considering the dynamical critical exponent Zi^t for integrated autocorrelations 
rather than for the exponential autocorrelation time; this is reasonable because it is the 
integrated autocorrelation which determines the cost of a Monte Carlo computation. 

1 Introduction 



Stochastic overrelaxation, and especially its variant usually referred to as "hybrid overre- 
laxation," is generally considered to be the most efficient algorithm for generating lattice 
configurations in pure gauge theory. As such it is also frequently used for the purposes 
of quenched simulations, although at the currently studied (small) correlation lengths the 
relative improvement over the standard Metropolis or heatbath local algorithms is a priori 
rather modest. 

The idea of generalizing overrelaxation methods to the stochastic case is due to Adler [2], 
who proposed an algorithm for multiquadratic actions such as free field theory (see also [3]). 
The initial expectations of improved performance were confirmed by studying the dynamical 
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critical behaviour of the algorithm in this context analytically [4-7] . The main result of these 
studies was that the algorithm can be tuned to achieve the dynamical critical exponent 

^cxp = 1 as compared to the generic value z^xp = 2 for the standard local Metropolis or 
heatbath algorithms. Some additional theoretical insights were also obtained in references [8- 
9]. 

Several extensions of the algorithm to interacting field theories and subsequent partial 
modifications were introduced by various authors, including [5, 10-14]. These algorithms 
were studied in more detail in several numerical works such as [15-18]. While these al- 
most invariably claim useful improvement, the truly systematic study (such that could make 
reasonable conclusions about dynamical critical exponents) is still missing. 

The purpose of this paper is to analyze the free field behaviour of overrelaxation further. 
Apart from our original aim of evaluating certain autocorrelations, we obtained some inter- 
esting conceptual insights which we feel are worth communicating. It has become clear over 
the years that the performance of Monte Carlo algorithm should be assessed with respect 
to a given operator fl (or set of operators): in particular, the most relevant characteristic 
then is the dynamical critical exponent zq, corresponding to the integrated autocorrelation 
for that quantity. While the dynamical critical exponent for the exponential autocorrelation 
time (which we denoted by z^^^ in previous paragraphs) usually gives us an upper bound, it 
often happens that zn < z^^p for most observables. 

If we have a situation where for operators of interest, za can be made smaller than 
z^^p, then the latter is somewhat irrelevant for the task at hand. In lattice field theory 
we actually have a specific group of operators we want to speed up, namely the ones for 
which the low energy (momentum) states are most important, because these operators are 
relevant for the continuum limit. We should thus attempt to optimize the overrelaxation 
algorithm for these quantities, rather than minimize ^^^^p. It has been implicitly assumed 
in previous studies that zero momentum operators (such as the magnetic susceptibility) 
actually saturate the Zexp-bound and so the above distinction is just a wishful thinking. Our 
main point in this paper is that while this is true for even/odd updates, where the dynamics 
of overrelaxation intrinsically couples the low and high frequency modes, it is not true for 
lexicographic updates, where the modes decouple at large volumes. In fact in the latter 
case we can tune the overrelaxation algorithm such that critical slowing down is completely 
eliminated for quantities that depend only on the zero momentum mode. We emphasize that 
this is a different tuning than the one leading to z^xp = 1; in fact, with this tuning z^^^ = 2. 
However, in this case the 2;exp-bound is actually saturated by high frequency quantities most 
of which we are not interested in. 

In the light of the above discussion, the "conventional wisdom" about the inability of local 
algorithms to perform better than z^^p ~ 1 (see for example [9]) can be rather misleading. 
While the statement as such might be true (but has not been proved even for free field 
theory), it does not tell us much about how the local algorithm will perform in the specific 
case of interest. 

The second point we want to make concerns the introduction of the additional noise in 
the overrelaxation update. In particular we attempted the same optimization for a scheme 
in which the updated site is chosen at random from a uniform distribution. It turns out 
that in this case it is impossible to tunc the overrelaxation parameter to reduce the critical 
slowing down for zero momentum quantities; in fact we cannot do better than z^ = '2 for 
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these operators. This means that overrelaxation can easily lose its magic if we introduce 
extra noise in the procedure. 

In [1] it was shown that the overrelaxation algorithm for free field theory is a special case 
of the class of Hybrid Molecular Dynamics algorithms where the fields are updated one site 
at a time using an analytic solution of the equations of motion. This rather suprising and 
elegant connection is especially intriguing in the light of the fact that, as shown in [14], a 
corresponding algorithm can be found for gauge theories. Here we adopt this point of view 
and will frequently refer to free field overrelaxation as Local Hybrid Monte Carlo (LHMC). 

We start in section 2 by formulating the LHMC algorithm and developing the techniques 
to calculate the integrated autocorrelations for typical zero momentum quantities, namely 
magnetization and susceptibility. In section 3 we analyze the performance of LHMC with 
even/odd, lexicographic, and random updates. We also give a proof that the integrated 
autocorrelation for the magnetic susceptibility vanishes at large volume with optimally tuned 
overrelaxation parameter when using the modified lexicographic update. In Appendix A we 
give an explicit asymptotic calculation of the update matrix for the lexicographic scheme, 
and in Appendix B we develop a formalism for bounding the finite volume corrections using 
methods of functional analysis. 

2 Local Hybrid Monte Carlo for Free Field Theory 
2.1 Analytic Solution 

We wish to consider a real scalar free field theory described by a functional integral with the 
usual action 

and a fiat (Lebesgue) measure for the 4> field. The lattice Laplacian operator is defined as 

d 

where /t is a unit vector in the /j, direction. Since the theory is free it can be diagonalized 
using the Fourier transformed fields 

where 0* = 0_p since (p^ G M.^ The action simplifies to 

s[4>] = I E fM' (1) 

^For notational convenience we will frequently omit the tildes and follow the convention that Fourier 
components are implied by subscripts p,q,... as opposed to x,y, 
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with 



cos ■ 



In particular we note that the lowest and highest frequencies of the system are /(o,...,o) — 
and ff-i r 1 r^ = + 4d. 



2.2 Local Hybrid Monte Carlo Updates 

Consider the dependence of the action S upon one degree of freedom only, S{(j).j_ 
^x<Px + constant, where here u;^ = + 2d and^ JF^. = l]^=i(0x+/i + (px-fi)- Introducing the 
fictitious Hamiltonian H{4>x,t^x) = ^t^I + S{4>x) we have 



H 



2 I 2 



+ constant, 



and the solution of the corresponding equations of motion is 

^x / , ^x\ J. , 7''a;(0) . 



(PxV) 5" = 0a;(O) 5- cosa;^ H — sma;i. 



or 



0; = (1 - 00. + 



C(2-C) , Tx 



c 



UJ UJ 

in terms of the overrelaxation parameter ^ = 1 — cosa;^. 

Assembling the field variables and their conjugate fictitious momenta into the vectors 



(2) 



/ 'Pxi \ 



= 



/ ^xr \ 



and TT = 



we can write the elementary local update (2) in the form 

^ M^cf) + P^tt; 



(3) 



(4) 



here Mj, and are x L*^ matrices, and the subscript refers to the site being updated. 
More explicitly, we have 



^x,y ~l~ ^x,z 



c 



UJ 



(5) 



^x,y^x,z 



a; 



^It is easy to get a factor of two wrong here! 
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Definition (3) implicitly assumes that we have ordered the field variables in a certain way. 
Sweeping through the lattice in this order, the complete update is given by 




E n 1 p^. 

i=l \j=i+l 



71 = Mcf) + Ptt. 



(6) 



Note that the form of the matrices M and P depends upon the predetermined order in 
which we chose to sweep through the lattice. In fact, as we shall discuss, the efficiency of 
the algorithm (optimized for a given operator) can depend crucially on this order. 

In addition to such fixed-order updates we will also analyze the algorithm in which the 
updated sites are chosen at random. The formulae (4-6) carry over unchanged to this case 
except that the updated site becomes an uniformly distributed random variable over Z^. 



2.3 Cost of Measuring the Operators 

If we are interested in measuring the expectation value of some operator Q, and is the 
average over a sequence of T configurations generated by some Markov process, we expect 



Here Varl] = {Vl?) — is the intrinsic variance of O, and Aq is the integrated autocorre- 
lation defined as 



where is the connected autocorrelation function. The cost of measuring (VL) to a given 
accuracy depends upon Aq,, and for local algorithms (such as LHMC) this is the only rele- 
vant characteristic which we need to compute in order to ascertain the performance of the 
algorithm. For a large system near criticality 

^n~r" (l«e«i>), 

where is the dynamical critical exponent, corresponding to quantity O.'^ 

In this paper we are interested in studying the performance of LHMC when applied to 
measuring some interesting operators in the context of free field theory. To this end we shall 
mainly consider the magnetization M. and magentization squared Al^, where we define the 
magnetization as 



^Note that our definition of the dynamical critical exponent assumes that the large volume limit is taken 
first. 
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The expectation values of powers of M. are given by 

= = -yx = ^(^) ( - it „ I z 

J_ooC^0oe ^/o-Po [ r(i) i/jj It 2^ e ^• 

In particular, {M) = 0, (A^^) = ^ and (A^^) = We thus expect that the measured 
values will satisfy AT = ± i-^MM+i, and AP = ± ^ V^^^^- 



2.4 Autocorrelations 

The general formalism for calculating the autocorrelations of any operator that is polynomial 
in the fields has been developed in [19]. Here we shall use the main ideas of that approach 
and derive explicit formulae for operators we are interested in. 

For convenience we will work in momentum space. In particular, we assume that the 
momentum-space representation of matrices M and P that characterize the updating scheme 
of our choice is known (see equation (6)) and we express the autocorrelations in terms of the 
corresponding matrix elements. 



2.4.1 Magnetization 

If we are interested in the Monte Carlo evolution of quantities linear in field variables, we 
can average the relation (6) over the independent fictitious momenta, and write 

Consequently we have 

l + AM = tc'M{t)-t ^^^'^^2^^'^ = = (I - M)o- , (7) 

t=0 t=0 vPol t=0 

where we have used the relation 

kMi) = —75-- 

An explicit representation for the normalized connected autocorrelation function (t) 
may be obtained by introducing an auxilliary variable x and evaluating the generating func- 
tion 

00 00 
1 + Am{x) = Y.^'CU{t) = E(^^)o,o = (I - xM)-,l. 

t=0 t=0 

On the right hand side we have a rational expression Q{x)/ R{x), where the degrees of 
polynomials satisfy the conditions^ 

deg R < L^, deg Q < deg R. 



the matrix (I —xM) ^ is block diagonal then degi? is equal to the dimension of the block containing 

the element (I — xM)^^. 
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We can thus perform the partial fraction expansion of Q{x) / R{x) and write^ 

oo degR degR oo 

t=0 j=l ^ ^3-^ j=l t=0 

Equating the coefficients of we can write the autocorrelation function in terms of the 
coefficients in the partial fraction expansion, namely 

dcg R 

Note that this also allows us to express the cumulative autocorrelation function in the form 

T degi? 1 l(T+1) 

t=o j=i 
and the integrated autocorrelation as 

2.4.2 Quadratic operators 

Turning now to quadratic operators we consider a generic quantity of the form 

n^^cj>;ci>,K, (8) 

q 

with the spectral density Kg being some function of the momentum. The change of the 
quadratic monomial (j)*4>q after a single sweep is given by 

r,s 

which can be averaged over fictitious momenta to give 

r,s r 

We can express this as 

= + (9) 

where 

<t>lq ^ <t>;<t>q. Mg^^^^ ^ (M* M),q,rs = M;,M,,„ and P^«, ^ E ^pV^,- 

r 

^The generalization to the case where R has degenerate roots is trivial. 
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The integrated autocorrelation is then given by 



" ho Egi, 92^91 ^92 [(^^l 051 092 052 )- (0^1051) (0^2 052)] 

Performing the sum over t and using the relation 

{4>*M*r4>s) - (0;05)(0:0.) = '^KA^sf,^ 

we find 

Note that the matrix elements of the identity matrix in the quadratic basis are Ipq^rs = ^p,r^q,s- 
In the case of M"^ we have Kg = S^^q, which imphes 

l + ^M2-[I-M%-o%. 

Note the formal similarity of this expression to that of equation (7) for the magnetization. 

In what follows we will also refer to the integrated autocorrelation for the energy (ac- 
tion) E; in this case we have Kg — giving 

^ + ^E=^,U^- M%lgg. 

There is obviously a certain redundancy in the matrix M^, since the basis elements we 
have used are not independent. In particular 0*0q = {4>*4>p)*. In practical calculations it 
is usually advantageous to reduce the basis to its independent subset, thus reducing the 
dimensionahty of the matrices involved. Also, in case of A4'^ we can use the basis {0p0g} 
instead of {0p0g}- 

Notice also that the formalism for calculating autocorrelation functions which we de- 
scribed in detail for the magentization carries over unchanged; calculating the involves 
the partial fraction expansion of ^^2(2;). 



3 Three Updating Schemes 

The exact analysis of autocorrelations for the linear stochastic update (6) is rather intractable 
in the general case, i.e., for an arbitrary ordering of the updated sites. However, in special 
cases the equations simplify and an exact analysis can be performed. In what follows, we 
will analyze the even/odd, lexicographic, and random updates. 

Prom our discussion in previous sections it follows that our work will split into two 
steps: Pirst we need to find the Pourier representation of matrices M and P corresponding 
to the updating scheme in question, and then we must evaluate the relevant formulae for 
autocorrelations. Note that the explicit form of matrix P is not needed for quantities we are 
interested in; this, however, is not true in general. 
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3.1 Even/Odd Updates 

Consider splitting the lattice into odd and even sites, as on a checkerboard. If we choose 
to update all the sites of a given sublattice before the sites of the other sublattice, then the 
result docs not depend on the order of sites within the sublattices. Indeed, updating all the 
even (odd) sites we may write 



(t>'x^(t>x+ X±{x) { -C(t>x + ^ —T^x + ) , 



1 ± (-1)^^=1'^'' 

The function x±(a;) has a simp' 



where x±(^) 



which is one for even (odd) sites and is zero otherwise, 
e Fourier representation. 



e I- ± e i 



where Pc = i(-^, • • • , -^), so we obtain in Fourier space^ 
0p = 



( Jc{2 - C) 

(/•p + 2 j ~C(0p ± 0p-Pc) H (t^p ± TTp-pJ 

H — o 2^ cos — T cos 



In matrix notation we have 



^(^P ^ + f^P J [ <PP-Pc ) 



C (2-C) / 1 ±1 
2a; V ±1 1 



7r„ 



P-Pc 



where /3± = (-1 ± ^ EJ=i cos ^) . 

For an even update followed by an odd update, which is the fundamental ergodic Markov 



^Note that cos M£z££k = _ cos 



9 



with 



and 



step, we have^ 

P^iPp - p;) 1 + Ppi^ + P; - p;) ) ' 

_ VC(2-C) /2 + /3+-/3- \ 

^" 2u; V Pp-Pp '^ + Pp-P^)' 

Note that with even/odd updates the Fourier modes of the field are only coupled in pairs. 
In particular, the lowest energy mode 0o is only coupled to the highest energy mode (pp^; for 
this case we will abbreviate Pq = — ^C{~^ i 2d/a;^). 

3.1.1 Autocorrelations for A4 

Using the formula (7) we can now trivially calculate the integrated autocorrelation for A4, 
namely 

Since ( G [0,2], this attains its minimum (as a function of () for ( = 2, where Am = ~h 
giving = ± 0. This is of course true even though, or rather because, C = 2 does not 
correspond to an ergodic algorithm. If we tune the overrelaxation parameter to ( = = 
(m^ + 2d)/{m^ + d) we have Am = 0, implying Zm = while the algorithm is ergodic. 
It should be emphasized though that this does not necessarily mean that the algorithm 
generates an independent estimate of A4 after every sweep: indeed, let us calculate the 
autocorrelation function According to our general discussion of autocorrelations, this 

requires the partial fraction decomposition of 

(T[-^M - 9M. - '^O + _ % . 



^If we are careful we observe that the equations of motion for the fictitious momenta are 
< = TT. + x± (x) I -Ctt. + ^^^^ ~ (-a;^. + J^) I 



or 



but 



) ^ ( vvf ) ( ) 4 1 ) ( .1 

Ti 1 ; V ^p-po ) V =Fi 1 y V ^p-p. 

so it makes no difference whether we evolve the momenta or not. Of course, this just reflects the fact that 
the momenta on even and odd sites are independent random variables. 
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in particular, we have 



Qo = i, gi = -(r + 

Ro = 1, i?i = -2(/3+ + /3- + 1) - (/?+ - 
If Xi and a;2 are tlie roots of the equation R{x) = 0, then 
Qo + QiXi Qo + Q1X2 



i?2 = + /3- + 1)^ 



R2{X2 — XijXi 



02 



i?2(a;i - X2)X2' 



Xi 



1 

X2 



(12) 



(13) 



Note that since X2 = a;*, we have a2 = a* and 62 = 61 as one would expect since the 
autocorrelation function Cj^{t) = aib\ + 02^2 is real. 
For ^ = we have in particular 



m 



02 = a^, 



(-\/d — i 



im 



|2' 



thus the autocorrelation function has an oscillatory behaviour with regions of correlation 
and anticorrelation. It is only the integrated autocorrelation that sums to zero for any m, 
giving zm = ^- 

3.1.2 Autocorrelations for M'^ 

The update matrix M is block-diagonal, so it suffices to consider the 4x4 block spanned 
by the basis {4>p4>p, 4>*(j)p-p^, (j)*_p^(j)p, (j)*_p^(f)p-p^} . The matrix elements of M and the above 
basis are real so this can be further reduced, and we use 



p 



where we have abbreviated {Mp)ij = Mij] after a straightforward calculation we obtain 





(l>p(pp ^ 




( Mi, 


2M00M01 


Mi, 




p^p-Pc 




MooMio 


MooMn + MoiMio 


MoiMn 




-pc^P-Pc / 




^ M?o 


2M10M11 


Mh 



(I - M«)oi 



1 



+ 



2/3; - - 2 



4(/3+ + /?- + 2) ' 8/3+ 4(/?+ + /?-)[(/?+ + 1)2 + (/?- + 1)2]- 

The integrated autocorrelation Am'^ — Mq)q I — 1 is then given by 

/ (m^ + Adm'^ + 9d^m'^ + 4rf3)(4 
-2(3m6 + IScim^ + 23d^m'^ + 12d:^)C 
+ {13m^ + Q2dm^ + 102^2^2 + 56^3)^2 
-4(3m6 + IQdm^ + 28^2^2 + lQd^)C 
+4(m^ + Qdm^ + 12ci2m2 + 8^^) 



(2 - OC^' 



/ (m^ + 4rfm2 + 8d2)(2 

-4(m*^ + Adm? + 4(i2)( 
V +4(m^ + 4dm2 + 4(^2) 

This equation indicates that without tuning of the overrelaxation parameter wc have z^2 = 2 
as expected. To find an optimal tuning we need to minimize Am'^ with respect to (. Since 
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we are only interested in the asymptotic behaviour as m — > — > oo), it suffices to consider 
just the leading terms as this limit is approached; in particular, we have 

d{2-C) , 7C^-16C + 8 , 



which is minimized by 



8C(2 - C) 



m 



2--^ + 0(m'). 



With this tuning the integrated autocorrelation becomes 

Vd 
2m 

giving zj^2 = 1. It is not possible to achieve 2; = as in the case of magnetization. 

For completeness we give the expressions for the autocorrelation function. Following our 
general strategy we find 



Qo + Qix 



^ ° ^^'^ " t1 1 - ~ Ro + Rix + R2X^ 1 - 63a; ' 



with 



Qo = (/3+ + 2)2 + /3-(/3- + 4), Qi 



+4(3-^(3+^ 

-(i + 2r'(i + 3r)(i+r))/5+' 
-(r + ir(r+2r 



and 



Ro = (/?+-r)'+4(i+/?++r), 



Ri = 



+4(1-/3-)/?+' 
+2(3 - 2/3- + 3/3-^)/3^ 
+4(1 -p-- p-' - /3-')/3^ 



V +(2 + 4/3- + 6/3-' + 4/3-' + /3-^) J 
R2 = i?o(/3+ + /3- + l)^ 

The coefficients ai,6i,a2,&2 are obtained by inserting these expressions into equation (13), 
while 03, 63 are given by 

2/3+/3- , 



3-3 



03 



(l + /3+ + /3-)^ 



In Figure 1 we plot the autocorrelation function and the cumulative autocorrelation 
function for Al^ in two dimensional free field theory, together with numerical data for various 
sample sizes. The error bars on the numerical data were obtained by binning n measurements 
into ^/n bins. The difficulty in estimating the integrated autocorrelation from small data 
sets is immediately apparent. 
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Figure 1: Comparison of theoretical prediction with measured values of Am^ two- 
dimensional free field theory using even/odd updates. The mass m — 0.2 and the over- 
relaxation parameter ^ = 0.95. 

3.2 Lexicographic Updates 

In this section we analyze updates in which the variable (px-fi is always updated before the 
variable (f)^. Strictly speaking such an updating scheme does not exist on a finite lattice with 
periodic boundary conditions; for example, if wc start with variable 0i on a one- dimensional 
lattice with L sites, then the variable = (f)L will certainly be updated after it. Allowing 
for these violations on the boundaries results in the class of updates that we call lexicographic 
after its most common implementation. 

To find the update matrix in this case it is useful to write the linear stochastic update 
(6) in a different form. The local update (2) involves the updated variable and its nearest 
neighbours; some of the neighbours, however, might have already been updated in the current 
sweep. We may separate the dependence on the "old" and "new" variables and write 

0' = 00 + Ncj)' + Ptt; 
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note that the matrix P — \^^JQ{2^^C) is different from the matrix P of equation (6), in fact 
we have 

M = {l-N)-^0, P={l-N)-^P. (14) 

It is straightforward to find an exphcit form of the matrices O and N for the lexicographic 
ordering. We have 

Ox,y — (1 ~ C)^x,y H 2 53 ^x,y-ji + 'T^^'V 
^ d ^ 

where we have separated the translationally noninvariant part 

= 2 "^X^x^y+ji^x,,,,! ~ Sx,y-p.Sxij.,L)- 

The matrix R represents the violation of our ordering rule on the lattice boundaries; as such 
it constitutes only a correction to the translationally invariant part. 

To see this more explicitly and to take advantage of the translational invariance in the 
infinite volume limit, it is once again convenient to work in momentum space. We have 

O^O^ + yR, N^N^ - ]rR, (15) 

iJ iJ 

where the translationally invariant parts are given by 



and 



Rp,i = 4 E K,, (e-^'^^^'-/^ - e^-^"/^) ; (17) 



where in the last equation we have defined 5^^ = Y\,y^f^ ^p^,qi.- momentum space the matrix 
elements of the translationally noninvariant contribution are explicitly suppressed by a factor 
of 1/L: this is expected since the violations have support only on the boundary whose size 
relative to the bulk volume is given by this factor. 

A similar suppression should also be explicit in the form of the update matrix: in Ap- 
pendices A and B we show that 

M^|:iG<0^o(i(ij)]. (18) 

Here C/^^ < ^/d for m > 0, so that the neglected correction is exponentially small in lattice 
size. The matrix elements of G'-'-' have an L-independent bound and are of the form 

{m1v,w} 
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with the modified (^-function is defined by 

C'-'"'^ n ^p... (19) 

where {/^i, . . . , /x;} is a subset of the integers {1, . . . ,d}. The leading term of equation (18) 
is translation invariant and corresponds to the situation with no violations of the ordering 
rule, namely 

- < - [(I - = ^v. ■ (20) 

The 1/L corrections arise from ordering violations on sites which have exactly one component 
having the boundary value (1 or L). We still have translational invariance in the remaining 
directions, which is expressed by the presence of 5^ explicitly we have 

= E ^'vA (-C + ^ E cos , (21) 

with defined in equation (43) of Appendix A. Similarly the 1/L' corrections originate 
with the sites with I components taking boundary values, and the product of 5-functions 
ensures the residual translation invariance va. d — l directions. 

3.2.1 Am. and Am'^ ^ Infinite Volume 

The analysis of autocorrelations becomes straightforward in an infinite volume because of 
the diagonal nature of the update matrix. For the magnetization we have 

, ,x 1 1 + 2d d , . 

l + ^^ = (I-M).-,; = ^-^ = ^-^--,. (22) 

Note that this result is identical to that of equation (11) for even/odd updates, and choosing 

we have again Am = 0, giving zm = 0. Nevertheless there is an interesting difference: the 
autocorrelation function takes the form of a single exponential, namely 



CMit) = {M^,or = exp (tln \^+Jjf ) , 



giving the exponential autocorrelation time 
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tm vanishes at C = C^, consequently each sweep generates a new configuration in which the 
value of magnetization is completely decorrelated from the previous value. 
For we have similarly 

1 + Am^-{1-M )oo,oo - r^lM^2 - C(2-C)mV^ + 2d)' ^^^^ 

and obviously TA42 = tm/'^- The integrated autocorrelation has a minimum at Q as defined 
in (23), where Am'^ vanishes, = 0, and complete decorrelation is achieved. 

We have thus obtained the interesting result that with a lexicographic ordering of sites it 
is possible to tune the overrelaxation parameter so that critical slowing down is completely 
eliminated. This is true for any quantity f2(0o) which depends only on the zero momentum 
component of the field: the Monte Carlo evolution of 0o is not coupled to the other modes 
in the infinite volume limit, and for one sweep we have 

(t>o^<t>'o = M^^o<Po + P^o^o, (25) 

where = ^^C(2 - C)(I - N^)"^ (from equation (14)). Since M^q = at C = Cc we get 

(n((/)o)n(0')) = (O(0o)^^(Po>o)) = (O(0o))(^^(0')), (26) 

thus the autocorrelation function as well as the integrated autocorrelation vanish. 

This is in contrast to the situation for even/odd updates: for lexicographic updates the 
Fourier modes decouple completely (in the infinite volume limit) and each mode is updated 
independently. Consequently, the autocorrelations of any quantity which depends on a single 
Fourier component can be tuned in an optimal fashion. For the even/odd case the modes 
are coupled in pairs and critical slowing down for, e.g., can only be reduced to zm'^ — 1 
at best. 



3.2.2 Other Quantities at Infinite Volume 

Let us now try to understand in a little more detail what happens at C = Cc- Consider the 
"staggered magnetization" squared, 0^^, which is a function of only the highest frequency 
mode, at C = Cc we have 

_ (m^ + 2dr 
^•^'^-~m2(m2 + 4d)' 

and thus z^2^ =2. In other words at C = Cc overrelaxation does not help if we want to 
measure this is not surprising as we have not sped up the high frequency modes; we 
took advantage of the fact that lexicographic updates allowed us to ignore these modes in the 
situation where we were not interested in them. On the other hand, if we want to measure 
(pp^ we can tune C to Cc = {m'^ + 2d) / {m^ + 3d) and eliminate critical slowing down completely 
for this quantity. The even/odd updating scheme does not have this fiexibihty. 

If we are interested in measuring both 0q and near criticality we are better off per- 
forming two lexicographic simulations without critical slowing down than one even/odd 
simulation. This is an extreme case of what we believe is a generic feature: if we want to 
optimize simulations with overrelaxation updates for some set of operators we can usually 
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find a better solution with lexicographic ordering of sites than with even/odd ordering at 
large volumes. The reason is that the cvcn/odd ordering introduces more constraints into 
the optimization than lexicographic update by coupling the modes in pairs. 

It is not easy to find an interesting operator that cannot be optimized well with lexi- 
cographic overrelaxation. For the sake of simplicity let us stay in the domain of quadratic 
operators of the type given in equation (8). As we have already observed, if K{q) is peaked 
at any single momentum wc expect ~ 0. If the operator Q couples to both low and high 
frequency modes, the contribution to autocorrelations of high frequency modes relative to 
low frequency ones will be suppressed by the fq factors in equation (10). To illustrate this, 
consider the simple quantity fl = (1)1 + 0^^; using formula (10) we obtain 

1 1 m'^ 

+ 



1 , ^^^^ _ 1 - «o)^ 1 - i^' + _ 

^ (m2 + 4dy 

We obviously get = at C = (c because the the 1/m^ behaviour of 1/[1 — {Mp^^^J'^] there 
is cancelled by the factor from its coefficient. Examples of important quantities of this 
type are two point functions. 

According to the above argument, the most "dangerous" quantities are the ones for which 
K{q) favours high frequency modes but does not completely decouple from low frequency 
modes. The obvious such operator is the energy, where the spectral weights exactly cancel 
(see equation (1)). In one dimension we get ze = 1 without tuning, ze = 2 a.t ( = (c, and 
C can be tuned to get ze — 1/2 for C = 2 — 2-\/2m; optimization for even/odd updates also 
leads to 2;^ = 1/2. 

Note that according to the above discussion, we should get Zj^t ~ for most of the 
quantities of interest for ^ = which optimizes functions of the zero momentum Fourier 
component. Furthermore, if we insist on measuring many operators with difi'erent spectra 
in one simulation, we can interleave lexicographic sweeps with C — Cc with sweeps at ( 
corresponding to z^xp = 1. This will give Zi^t — for most of the important operators while 
ensuring that Zi^^ < 1 for all operators. 

3.2.3 Autocorrelations in a Finite Volume 

While the value of 2;;,^ is defined by the infinite volume behaviour, it is useful to understand 
the finite size corrections to autocorrelations. Interestingly, for the magnetization there are 
no corrections to the infinite volume result: indeed, we have 

(I - M)-^ = [l - (I - N)-^0]~^ = (I - AT - 0)-\l - N). (29) 

From equation (15) we have I-A^-C> = I-A^^- O^, so 

(I - M)-i = (I - M^)-^ + y{I -N"" - O'^y^R] 

since i2o,o = 0, it follows that 

1 + ^A, = (I - Uyl = (I - M^)o-J, 
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and the result (22) is exact for any volume. 

In case of Ai"^ the situation is more complicated. (I — M'^)~^ is expected to have a 
1/L correction to the translationally invariant part. We have not been able to evaluate the 
correction in the closed form and we have no reason to believe that the zero momentum 
matrix element vanishes; in fact, our numerical experiments in one dimension confirm that 
the correction is of order 1/L as is illustrated in Figure 2. The data shows a very small 
curvature at large values of L, which corresponds to the small negative intercepts for the 
linear fits; we do not understand the cause of this subtle effect. 




0.00 0.02 0.04 0.06 0.08 0.10 0.12 

1/L 



Figure 2: Am^ for lexicographic updates in one dimension with m — 0.125 and 0.250 at 
C — Cc- The error bars are smaller than the symbols. 

3.2.4 Random Lexicographic Updates 

On physical grounds it is obvious that the effect of surface terms on behaviour of An for 
any quantity Q will go to zero, typically as 1/L, in the limit of large volumes. A periodic 
lattice does not have any boundaries, it is only the updating scheme that introduces a 
boundary as the collection of points at which the lexicographic property is violated; for this 
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reason expressions describing autocorrelations do not have exact translational invariance. An 
obvious way to enhance translational invariance is to randomize the choice of starting point 
for each sweep through the lattice. For this case we are able to provide a relatively simple 
direct proof that surface effects only give small corrections to the integrated autocorrelation 
even for quadratic operators. 

The boundary is completely specified by the site that is updated first in the sweep; if 
X = {xi,X2, ■ ■ ■ ,X(i) is this first site then a site y belongs to the boundary if i/j — Xj or 
Vj ~ ^j+L-i foi' at least one j. We thus need to give every lattice site an equal chance of 
being the starting point of the lexicographic chain. 

Consider the representation of lexicographic updates in terms of the local update matri- 
ces (6) 

where the sequence Xi, 2:2, . . . , x^a defines the lexicographic ordering. If we relabel the sites 
by a lattice translation x ^ x — a and apply the same lexicographic update in the new 
coordinates then in the old coordinates this corresponds to choosing the initial site to be 
xi + a and the update matrix to be 

The matrix M^+a can be written in the form 

Using the composition property S'^"'^S'^^'> — S^"''^^'> we have 

Mia) ^ 



We now propose the updating scheme where before any given sweep, this relabelling 
of sites is performed at random. In other words, the translation vector a will be a ran- 
dom variable, uniformly distributed over Z^. The linear update matrix averaged over these 
translations is given by 

M = ^ y S^'''>MS^-''\ 

Td ^ 

In Fourier space we have 

5'(«) = S e~'^P'"' and M ^ S M ■ 

as expected, the linear update matrix retains only the diagonal part of M. 
The quadratic update matrix in the {(t>p(t>q\ basis is similarly given by 

Performing the sum over translations yields 
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implying that the quadratic monomial 0*0^ couples only to monomials 0*0s satisfying p—gr = 
r — s. Restricting ourselves to the sector with p — q — and defining $p = (f)*(f)p we have 

$' = M«-o$, with M^f = M;^^M,,,; 

here M^'^ is the corresponding L'^ x L'^ block of the quadratic update matrix given by the 
Hadamard product M* ■ M. We are interested in showing that 

1 + Am^ = (I - M* • M)-J ^ (I - M^* ■ M^)-J as L ^ oo. (30) 

Let us start by writing the matrix M explicitly as (18) 

1=1 ^ 

where we have absorbed the exponentially small corrections into the matrices G^''^; we then 
have a decomposition 

I - M* • M = I - M^* ■ - E^^^ - E^^^ = I - M^* ■ - E, 

where 

1=1 ^ 1=1 ^ 

d -I d -| 

ir{2) _ n('-)*n(^) = V" sill,.. . M s:v]_,.. uj ritii,... ,1X1* f:<vi,...,uj (oq\ 

^p,q j^l^j^p,q ^p,q j^l^j^p,q "p,q ^ p,q ^ p,q ■ \^^) 

{ui,...,Uj} 

We now argue that the "error" matrix E is small at large volumes in the sense that its 
matrix norm tends to zero in that limit. To show this we use the maximum sum row norm 
\\A\\ = maXpX^q \^p,q\- Consider first E^^\ which is a diagonal matrix since is diagonal; 
the matrix elements of and G^'^ are bounded by L-independent constants, thus there 
is a constant Ci such that < Ci/L for all L. Turning to E^'^\ since there are an L- 

independent number of terms in equation (32), it suffices to consider the norm of a general 
term in that sum. We have 

j-l+j^p,q ^p,q ^p,q ^p,q ^ Jj+j-k — ^ ' ^ > 



q ^ 



where k is the cardinality of the set {/^i, . . . , /x/} fl {i/i, . . . , i/j}; we have again used the fact 
that the matrix elements of the matrices G/^i' - 'W have an L-independent bound iGp,^^ "''^' | < 
C. The last inequality in (33) follows from the fact that k < min(/,j) and consequently 
1 + i - k >l. Wc have thus established that \\E^'^^\\ < C2/L for all L. Putting the two 
bounds together we have 

" " - L ~ T' 
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As the final ingredient we use the standard formula for the error in the inverse, namely 



(A - E)-^ - A-^ 



< 



1 - \\A-^E\ 



\E\ 



(34) 



this is valid^ for any matrix A and "error" E and any norm in which < 1. 

Let A = I—M^*-M^, then \\A^^ \\ has an L-independent upper bound and limi^oo = 
0, so we have ||74~^£J|| < ■ \\E\\ < 1 for sufficiently large L. Formula (34) is thus 

applicable, and 



M)-^ - (I - ■ M 



(35) 



for sufficiently large L; the required result (30) follows. 

In the infinite volume limit A%i is the same as for the standard lexicographic update (24). 
Whilst equation (35) tells us that the leading finite volume correction is at worst 0(1/L), 
it is actually (9(1/L^) at C = Cc- Indeed, the matrix elements of E^'^^ are supressed by 
and consequently the 0{1/L) part is determined solely by E^^\ which is diagonal. Using the 
explicit form (21) of matrix G^^^ this leading correction can be evaluated, and the correction 
vanishes at C = Cc- The behaviour of the correction is in accord with numerical results, 
as shown in Figure 3. 



3.3 Random Updates 

Finally consider an updating scheme in which the sites to be updated are chosen at random. 
A sequence of L'' such local updates will be called a "sweep" . 

For a single site update we use the update matrix given in equation (5) with the 
updated site z being a random variable, uniformly distributed over Z^. In Fourier space we 
obtain 

W,,, = 5p,, + -^e-^-'(^-^W^/(g), 

with 



m = c 



2 ^ 27r 
— cos —q^ - 1 

n=i ^ 



For autocorrelations of quadratic operators we define the quadratic elementary update matrix 
in basis {4>p(l)q), namely 

iM%,,rs = M,)pg,rs = (M,);,(M,),,,. 

3.3.1 Autocorrelations for M. 

The linear update matrix for a random update sweep, M, averaged over the independent 
random variables Zi is given by 



^The derivation involves the expansion of {A — E) ^ in a Neumann series. 
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Figure 3: Am^ for random lexicographic updates in one dimension with m — 0.25 at C = Cc- 
The fit was made to the leftmost six points, but not constrained to go through the 

origin. 



Explicitly, we have 



giving the integrated autocorrelation 



1 C^m^ / 1 \ 

' l-exp(-JS,) L^8K + 2ci)^sinh^(,^) 

In the small mass limit we have 

2d 

indicating that zj^ = 2 regardless of the choice of (. Even though the modes have completely 
decoupled and the autocorrelation function is a single exponential, the random choice of the 
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site to be updated completely destroys the coherence and the algorithm cannot be tuned to 
reduce the critical slowing down. 



3.3.2 Autocorrelations for M'^ 

The quadratic update matrix M'^ averaged over the random variables Zi is 



^^ = n(i^E^S) = (i;^EMf) . (36) 



Straightforward evaluation gives 



^i:a4«=i+;^(c«+>«). 



with 

D^rs = ^AAfir) + f{s)), C^,,, = 5,.r,,-sfir)f{s) 
Inserting this into equation (36) yields 



M« = exp(L>«) + O (-^) ; 



we have not evaluated the 0{1/L^) correction in this case. In the infinite volume limit we 
obtain 



ipq,rs _ ^f{r)+f{s) ' 

and the integrated autocorrelation of M'^ is given by 

l + ^^. = ^---i-^«^ + 0(m°) form^O. 

We therefore find that again that the algorithm with random updates cannot be tuned to 
reduce critical slowing down below z_xi2 = 2. 

4 Conclusions 

The principal lesson which we may draw from this analysis is that it is possible to reduce the 
cost of a Local Hybrid Monte Carlo computation by a judicious choice of the order in which 
the variables are updated, as well as a careful tuning of the amount of randomness which is 
introduced into the system. While this differs from some previous claims, it is similar to the 
situation for deterministic Gauss-Seidel linear equation solvers. 

The relevant autocorrelation for the determination of the cost is the integrated autocor- 
relation, and it is in terms of this autocorrelation that we define the critical exponent Zi^t- 
Furthermore, we define z-.^t by ^ oc as L, ,^ ^ oo and ^/L ^ 0. This is to be contrasted 
with a finite size scaling analysis in which the limits are taken in a different order. For a 
zero temperature quantum field theory our choice of limits seems more natural. 
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Most local algorithms are special cases of LHMC, so our analysis is of fairly general 
applicability if we assume that the behaviour of autocorrelations in asymptotically free in- 
teracting theories is similar to that of free field theory. We intend to study this empirically 
in a future publication. 

Of course local algorithms are directly applicable only to theories with local bosonic 
actions, but they are used within some dynamical fermion algorithms such as that proposed 
by Liischer [20]. 

One issue of practical importance which has not been addressed at all in this paper is the 
difficulty of implementing the lexicographic scheme on parallel computers. Nevertheless it 
may well be advantageous to use a local lexicographic scheme on parallel computers, as has 
been suggested for the case of preconditioning conjugate gradient linear equation solvers [21]. 
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A Calculation of M in the Lexicographic Scheme 

Our aim in this section is to show that the update matrix M for lexicographic updates has 
the structure as given in equation (18), with the leading term of equation (20) and subleading 
term of equation (21). We shall prove this by induction in the number of dimensions d for 
fixed lattice size L. The update matrix is a function of the mass and the overrelaxation 
parameter as well as d and L, M^'^'> — M'^'^\m, Q. The induction step assumes that M^"^^ is 
known for all m > and < C < 2 and relates M('^+^)(m, C) to M'^'^\m', (). 

Consider the vector of field variables (p = (0i, 02, ■ ■ ■ , 0l) in + 1 dimensions, where 
we only keep the index corresponding to the new dimension d + 1, that is, (f)^ is itself a 
vector of field variables in d dimensions. We first find the matrix M^'^~^^\ corresponding to 
lexicographic update of the variables 0^ only. Separating the dependence on the "old" and 
"new" variables in this subspace as usual, we obtain 

0^ = o(<^v. + A^^'V; + i4(0.+i - 0.-i), 



or 

0; = M(<^)0, + (I - ArW)-i-^(0,+i - 0,_i). 

We emphasise that the system under consideration is 0? + 1 dimensional, and so a;^ = + 
2{d + 1); furthermore the update matrix M^'^'> for a d-dimensional subspace is the same as 
that for a d-dimensional system with mass m'^ — 2-|-m^. For convenience we define a variable 
t = C/^^; and we consider A^^"^) to be a function of t, and O^''^ and M^'^^ to be functions of 
both t and We can thus write 
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The full update matrix is the product (in lexicographic order) of the subspace update 
matrices. In particular, 



(d+i) 



where Sx^y = represents a translation in the new coordinate d + 1. 

In Fourier space it is convenient to introduce a = e^'^*/^, so we have 



,>{<i+i) 



M 



1 



p,q 



M^'^ -I + t{l- ArW)-l(Q;« ^ ^-q^J I 



Note that we only show the indices corresponding to the new dimension. Writing T = 
t{I — A^^'^))"^ one can show by induction that the matrix elements of the Ith power for / < L 
are given by 



with 



CAd+i) ' 



p,q I ' L J J 



Vl{xi, ...,Xn) 



E 



^mi r. 
•^1 ■ ■ ■ -^r. 



m\ ,...,mn'>0 
^iH \-mn=l 



which satisfies the following useful identities 



j=0 

Vl{xi,...,Xn) = XnVl-l{xi,...,Xn)+Vl{xi,...,Xn-l) 



'Pi{x,y) 



^l + l _ yl+1 



x-y 

The Lth power has one extra term giving 



M('^+i)1_ = Up^q + ^ bL-i(Qi^ T) + T] [m('^) - I + i(I - Ar('^))-i(Q;« + a-^)] . (37) 

" " ^ L J L J 



P,9 



We can obtain a simpler recursion relation for (I — A^) ^, which is related to M through 
multiplication by the diagonal matrix; indeed, using equations (29) and (15) we get from (37) 



(38) 



Finally we need to transform VL-i{oi^, ofl, T) into an explicit form; straightforward algebraic 
manipulation leads to 

VL-i{a^, a\ T) = 5p^qL{an -T)-^- T{an - T)-^{an - T)-^{1 - T^). (39) 

RecaUing the definition of T we have also 



{an - r)-i = - 



ta-P 



I 



1 - ta-P' 



(40) 
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It is important to note here that the inverse on the right of this equation exists: we know 
that (I — N^'^^)~^ exists for all |t| < 1/d, where t enters as a multiplicative factor in definition 
of N^'^^ (see equations (15) — (17)). However, in 1 dimensions t satisfies a more stringent 
constraint, namely t <l/{d+l), and consequently 



t 



implying the existence of (I 
equation (38) we get an explicit recursion relation 



l-ta 



- ta-P 



< J- 
- l-t 

^ for all p 



< 



1 



Inserting relations (39) and (40) into 



1 - ta-P 



I- 



+ 



(I 



{l-ta-P){l 



ta-1) 



1 - ta-P 
I- 



ta-P 



ta-1 



(41) 



In principle this allows us to calculate the update matrix in any dimension. Consider 
building a 1-dimensional update matrix from single site updates {d — 0); Obviously, A^^°^ = 
and the recursion relation gives the exact result 



P,<1 



1 - ta-P 



t 



a 



-p 



l-t^) 



{I - ta-P){l - ta-1) 



(42) 



Note that the above result has the expected structure and the leading term exactly corre- 
sponds to [I — N'^^^ with N^'^'> denoting a diagonal part of N^'^\ defined in (16). Using 
we can calculate A^*^^) etc. The only problem from the practical point of view is the 
calculation of = t^[(I - N'^'^^)-^]^ . However, since t < 1 this term is exponentially small 
in L and can be neglected at large volumes. Similarly, the exponentially small terms will be 
generated (and can be neglected) from matrix multiplications in the recursion relation when 
using, e.g., the Poisson resummation formula for the matrix elements of the product. 

While the above program can be carried out to arbitrary number of dimensions, it is 
sufficient for our purposes to deduce the structure that follows from equations (41) and (42). 
Indeed, one can easily verify the following conclusions: 

1. In any number of dimensions we have 



The matrix elements of K^^^ are bounded by an L-independent bound and have the 
form 

where the modified delta function is defined in equation (19), and {/Xi, . . . repre- 
sents any subset of integers {1, ...,(/} with I elements. This follows from the (42) and 
the fact that the above structure is invariant with respect to recursion relation (41). 
Note that the exponentially small terms could also be included in matrices K^^^ ^ but 
we keep them separately for the practical reasons mentioned above. 
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2. The invariant form of the leading term is given by 

1 



which corresponds to [I — N^]~^ as expected. 
3. The invariant form of the first finite volume correction is 



■ (43) 



The corresponding forms for the update matrix M are given in equations (18-21). 

B Boundedness of Finite Volume Corrections 

In this appendix we show that the surface corrections for the lexicographic updating scheme 
do indeed go to zero as the inverse of the lattice size. For the case of quadratic operators wc 
shall assume that the infinite volume limit of det(I — M* ® M) 7^ 0: we have not yet found 
a simple proof of this, although it would follow trivially from a proof of the ergodicity of the 
LHMC algorithm in an infinite volume. 

The basic idea of the proof is to expand about the infinite volume limit, and thus our proof 
is naturally expressed in the language of functional analysis. Our argument is essentially 
that the update operator M is a compact (completely continuous) operator on the Hilbert 
space -^^^([0, 27r]) — in fact in one dimension it is a Hilbert-Schmidt kernel — and that the 
finite volume results may be expanded about the infinite volume ones using the Poisson 
resummation formula. We have chosen to prove our results following the original method of 
Fredholm, as this seems to be the most direct approach [22]. 

We start by proving a simple bound on determinants: 

Lemma 1 (Hadamcird) For any matrix A : ll{C} — > ll(C} 

idet^i <niiAii 

i 

where Ai is a row of A and — J2j lAjP ^■s its norm. 

Proof. If A has any row which is zero then the inequality is trivially satisfied. Otherwise 
construct a new matrix B whose rows are Bi = Ai/\\Ai\\, then the bound is equivalent to 
I det5| < 1 for all complex matrices whose rows are normalized. If the rows of B are not 
all orthogonal then there are two rows u and v such that |(M,f)| > 0, = \\v\\ = 1. If we 
replace the row v by v' = v — {u, v)u then the determinant is unchanged, and {u, v') — 0. On 
the other hand = \\v' +{u,v)u\\'^ = + |(m, f )P||m||^, so = l — \{u,v)\'^ < 1; thus 
replacing the row v' by v" = v' /\\v'\\ increases the determinant by a factor of l/H^^'H > 1, and 
all the rows of the resulting matrix have unit length. This proves that the maximum value of 
the determinant must occur when the rows are orthonormal. If the rows of B are orthonormal 
then BB'f = 1 and hence det BB'' = det B det = det B{det B)* = | det = 1. ■ 
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Corollary 1 // all the elements of the matrix A are bounded, \Aij\ < C, then IdetAj < 

{cVkf. 

Theorem 1 (Fredholm) Let A : Z|,(C) — > ^|,(C) he a matrix all of whose elements are 
hounded by some (L -independent) constant |Ay| < C; then det ^1 + 
a constant which is independent of L. 



< C where C is 



Proof. Let us introduce a parameter A and then expand in powers of A to obtain 

/ A,; 



det f 1 + ^A 



using Hadamard's bound 



E det 

ii>i2>--->ifc 



Ann \ 



ikik I 



k=0 



k 



{cVkr < E 



(ACV/O 



fe=0 



A;! 



C'{X) 



where the ratio test may be applied to show that the series converges for all A. The desired 
result follows upon setting A = 1. ■ 

Definition 1 The minor A^*''^ of a matrix A is defined hy 

1 if k — i and j — I 



^kl 



Afei ifky^i and j ^ I 
otherwise 



Corollciry 2 (Fredholm) Let i^he a matrix as in Theorem 1, then 



det(l + -A(*^) 



< 



C ifi^j 
C'/L tfi^j 



where C is an (L -independent) constant. 

Proof, li i = j then the matrix I + j-A^'^^^ is of the form to which Theorem 1 is applicable. 
If i 7^ J then introducing the parameter A and expanding in powers of A as before we obtain 

^ Aji Aji^ . . . Ajjj, \ 



det ( 1 + ^A(*J) 



fc+i 



E 



-1)'-+^ det 



fe=0 



n>i2>"->jfc 

ir>i>ir+l 

is>j>is+i 



using Hadamard's bound 

L-2 / ^\ k+1 



E {C^/kTl: 



fe+1 



^ikik / 



Sir 



4l>i2>--->«fc 
ir>«>«r + l 

is>i>js+l 



fe=0 



^ ^ ^ (ACx/FTT)'=+^ _ C"(A) 
— T ^ 



A:=0 



k\ 



where the ratio test may be applied to show that the series converges for all A. The desired 
result follows upon setting A = 1. ■ 
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Definition 2 The adjoint Adj M of a matrix M is defined to be the transpose of the ma- 
trix whose elements are the determinants of the corresponding minors, that is (Adj M)ij = 
detM(^'). 

Theorem 2 Let A : l1{C) — > i|,(C) be a matrix all of whose elements are bounded by some 
constant |Ajj| < C, and det (^I + |;A^ > a > is bounded below by an (L -independent) 

constant, then {\-\- -^A^ =1+ j^G, where all the matrix elements of G are bounded by an 
(L-independent) constant, \Gij\ < C" . 

Proof. From the algebraic identity (Cramer's rule) 

(l + A^ Adj {l + -^A^ = det {l + -^A^ I = Dl 

and the fact that [Adj (l + ;iA)] = C^/L for i ^ j and [Adj (l + ^Ajj = Cjj, which 
follow from Corollary 2, we have 



(l + Ia,,) C,, + Y: ^,^J^C^J = D. 



From this we find that Cjj = D + 0(1/ L) ^ Cjj/ D = 1 + 0(1/L), where we have made use 
of the fact that D is bounded below by the constant a. The desired result then follows from 
substituting this relation into Cramer's rule above. ■ 

Lemma 2 If the matrix elements Ajj referred to above are taken to be the values A (^^, ^) 
of a Lebesgue integrable function then we can define the infinite volume limit of the Fredholm 
determinant, and if this limit is non-zero then there must exist a non-zero constant a which 
bounds it from below. 

Definition 3 Let a be a subset of the integers {!,..., d}, and S^^^ = O/x^a ^p^^q^^■ ^ matrix 
A has support on dimensions a if — A^q, with p and q being a d-tuple of indices, 
p= {pi,...,Pd). 

Notice that according to our definition A^'^^ is a diagonal matrix. 

Corollary 3 There is a permutation of the indices p for each matrix A^""^ such that it is 
block diagonal with #a x #a blocks. Each block satisfies the assumptions of Theorem 2, so 
the inverse may be expressed in block diagonal form where each block is of the form specified 
in Theorem 2. 

Lemma 3 If A^"^ has support on dimensions a and B^f^^ has support on dimensions (3, and 
all the elements of are bounded by an (L-independent) constant, \A'/^^\ < C and \B^^\ < C, 
then 

Y,L-*^A^;^L-*fBf = 

where has support on dimensions a\J (5 and \ W^^^'>\ < C . is the cardinality of 

the set a. 
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Proof. Consider the matrix product 

V B'^^^ = ^ a B Y\ 6 6 



The product of Kronecker tensors may be spht into four disjoint classes, a U /3, /3 — a, a — 
and an /3: 



n 1 

i/iGan/3 



It is immediately apparent that only for in the last set of indices (a Pi /3) does the sum over 
Qi^ have more than one non- vanishing term, and that for /i in the first set of indices (a U /5) 
the expression vanishes unless — r^. We thus have 

with \Wpr\ < C The result now follows from the observation that Li P — #o; + #/3 — 
^anp. m 

Theorem 3 Let A^°^ be a matrix with support on dimensions a all of whose elements are 
bounded by an (L-independent) constant, and furthermore 

lim det I I + V L-#"A(") | ^ 0; 

then ^ 

where G^'^^ its has support on dimensions a and all its elements bounded by a constant. 

Proof. We may order the subsets of {1, . . . , rf} lexicographically, that is a < /5 iff jj^a < 
or #q; = #/3 and minjgQ_^i < minjg^_Q,j. Let 7 be the smallest set for which A^'''^ 7^ 0. 
Then 

I + ^"*"A(") = (l + L-#^AW) 1 1 + (i + L-#^A(t))"^ Y i^"*"A(") I . 

The infinite volume limit of the determinant of this quantity is non zero, therefore the infinite 
volume limit of the determinant of each factor must be non zero. Prom Corollary 3 we have 
that (I + L-#tAW)-^ = I + L-#tgW with \Gpq\ < C, hence 

I+^L-#"A("M = 1 1+ (l + L-^^G^^)) ^L-#"A("M (l + L-^^G^^)) 

a^0 
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where we have used Lemmas 2 and 3, and the fact that Q;>7^Q;U7>7in the penultimate 
hne. The final line follows by induction on 7. ■ 

One application of these results uses the matrices A^, O, N^, O^, and R introduced in 
the body of the paper, and the matrix A defined as A = (I — N^)~^R. 

Corollciry 4 M has the form given in equation (18). 

Proof. Observe that 

I - AT = (^I - AT-D + }-R^ = (I _ AT-D) (^i + 1(1 _ N^'y^R^ = (I - AT^) (l + 1a) ; 

since there is a basis in which A^ is a strictly triangular matrix it follows that det(I — A^) = 
1 (VL), and thus we may conclude that limi_»oo det (l + ;i a) ^ 0. Therefore A satisfies the 

assumptions of Theorem 3, so — j^Aj has the desired form. Prom equation (29) we have 

M = I - (I - N)-\I -N^ -0^)^I - (1+ 1a^ ^ (I - Ar^)-^(I -N^ - O^), 

so we may apply Lemma 3 to show that M also has the desired form. ■ 

If the infinite volume Markov process has a unique fixed point then from equation (9) we 
obtain upon averaging over this fixed point distribution of 

(I-M«)(0«)^,, = P«, 

and the uniqueness of (0*')<^,7r imphes that limL^oo det (I — M* ® M) 7^ 0. As a consequence 
of Corollary 4 we may write 

I - M* ® M = (l - M^* (g) M^) (l + Y. L-*''A^°'^^ , 

where the matrix elements of A^"'^ are bounded by an L-independent constant, and thus from 
Theorem 3 and Lemma 3 we may conclude that a result analogous to Corollary 4 applies to 
quadratic operators. 
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